% ****************************************************************
% test association method

addpath ../
addpath ../export_fig/

load genotype.mat;
GENO = single(SNP');
clear SNP*;


% ****************************************************************
% generate expression data
[nsnp,nind] = size(GENO);
ngene = 200;
regulated_gene = 100;
cis_per_gene = 10;
causal_per_gene = 3;


CIS = logical(sparse(ceil((1:ngene*cis_per_gene)/cis_per_gene), ...
                     randsample(nsnp,ngene*cis_per_gene,true), ...
                     1, ngene, nsnp));

[gidx,sidx,~] = find(CIS);

num_causal_pairs = floor(causal_per_gene*regulated_gene);
causal_idx = randsample(numel(gidx), num_causal_pairs);
causal_pairs = [gidx(causal_idx), sidx(causal_idx)];

W_true = sparse(causal_pairs(:,1), causal_pairs(:,2), ...
           (-1).^(rand(num_causal_pairs,1)<.5), ...
           ngene, nsnp);

EXPR = zscore(single(W_true * double(GENO)) + 0.1*randn(ngene,nind,'single'),[],2);

% ****************************************************************
[Z,W,Wvar,InvTauSqW,Theta,ThetaVar,InvTauTheta,ELBO] = run_modQTL_association(EXPR,GENO,[],CIS,100,500);

% ****************************************************************
h1 = figure(1);
set(h1,'position',[.5 .5 600 600]);
subplot(2,1,1);
set(gca,'FontSize',20);
hist(Z,0:.05:1); axis tight;
title('all SNPs');
ylabel('Frequency');

subplot(2,1,2);
set(gca,'FontSize',20);
hist(Z(unique(causal_pairs(:,2))),0:.05:1); axis tight;
title('causal SNPs');
xlabel('Posterior Inclusion Probability');
ylabel('Frequency');
export_fig('sim09_fig01','-pdf','-a1','-m2','-transparent');


h2 = figure(2);
labels = zeros(nsnp,1);
labels(unique(causal_pairs(:,2))) = 1;
[fpr,tpr,~,auroc] = perfcurve(labels, Z, 1, 'xCrit', 'FPR', 'yCrit', 'TPR', 'NegClass', 0);
minval = min(min(fpr),min(tpr));

set(h2,'position',[2 2 400 300]);
set(gca,'FontSize',20);
semilogx(fpr, tpr, '-o');
xlabel('False Positive Rate');
ylabel('True Positive Rate');
title(sprintf('AUROC=%.2e',auroc));
axis([minval 1 minval 1]);

export_fig('sim09_fig02','-pdf','-a1','-m2','-transparent');



h3 = figure(3);
set(h3,'position',[3 3 400 300]);
lab = zeros(ngene,1);
lab(unique(causal_pairs(:,1))) = 1;
bh=boxplot(1./InvTauSqW, lab, 'notch', 'on');
for bi=1:size(bh,1),
    set(bh(bi,1),'linewidth',2);
end
set(gca,'FontSize',20);
xlabel('Causal gene'); ylabel('\tau_w^2');
export_fig('sim09_fig03','-pdf','-a1','-m2','-transparent');
